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Basis of the report 

a. With regard to the language, the international search was carried out on the basis of the international application in the 
language in which it was filed, unless otherwise indicated under this item. 

I I the international search was carried out on the basis of a translation of the international application furnished to this 
Authority (Rule 23.1(b)). 
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With regard to any nucleotide and/or amino acid sequence disclosed in the international application, the international search 

was carried out on the basis of the sequence listing : 
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filed together with the international application in computer readable form, 
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furnished subsequently to this Authority in computer readble form. 
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the statement that the subsequently furnished written sequence listing does not go beyond the disclosure in the 
international application as filed has been furnished. 

I I the statement that the information recorded in computer readable form is identical to the written sequence listing has been 
furnished 

I I Certain claims were found unsearchable (See Box I). 
I I Unity of invention is lacking (see Box II). 



4. With regard to the title, 

PC] the text is approved as submitted by the applicant. 

I I the text has been established by this Authority to read as follows: 



5. With regard to the abstract, 

I X I the text Is approved as submitted by the applicant. 

□ the text has been established, according to Rule 38.2(b), by this Authority as it appears in Box 111. The applicant may, 
within one month from the date of mailing of this international search report, submit comments to this Authority. 

6. The figure of the drawings to be published with the abstract is Figure No. 17 

I I as suggested by the applicant. Q None of the figures. 

pr| because the applicant failed to suggest a figure. 

I I because this figure better characterizes the invention. 



Form PCT/ISA/210 (first sheet) (July 1998) 



INTERNATIONAL SEARCH REPORT 



International Application No 



A. CLASSIFICATION OF SUBJECT IWIATTl 

IPC 7 GOlVll/00 



According lo International Patent Classification (IPC) or to bolh national classitication and IPC 



GB 01/03042 



B. FIELDS SEARCHED 



Minimum documentation searched (classification system followed by classification symbols) 

IPC 7 GOIV 



Documentation searched other than minimum documentation to the extent that such documents are included in the fields searched 



Electronic data base consulted during the intemational search (name of data base and. where practical, search terms used) 

INSPEC, COMPENDEX, EPO-Internal , WPI Data 



C. DOCUMENTS CONSIDERED TO BE RELEVANT 



Category * Citation of document, with indication, where appropriate, of the relevant passages 



Relevant to claim No. 



LEWIS ROLAND W ET AL: "FINITE ELEMENT 
ANALYSIS OF SURFACE SUBSIDENCE" 
INT CONF ON EVAL AND PREDICT OF 
SUBSIDENCE, PAP;PENSACOLA, FL, USA JAN 
15-20 1978, 
1978, pages 400-416, XP001036503 
1978 ASCE, New York, NY 
page 401 

US 5 417 103 A (HUNTER ROGER J ET AL) 
23 May 1995 (1995-05-23) 
column 6, line 32 -column 7, line 28 
column 13, line 52 -column 14, line 33 

-/-- 



1-4,7-10 



I, 5-7, 

II, 12 



Further documents are listed in the continuation of box C. 



Patent family members are listed in annex. 



* Special categories off cited documents : 

"A" document defining the general state of the art which is not 

considered to be of particular relevance 
*E' earlier document but published on or after the international 

filing date 

'L' document which may throw doubts on priority claim(s) or 
which Is cited to establish the publication date of another 
citation or other special reason (as specified) 

'O" document referring to an oral disclosure, use. exhibition or 
other means 

'P* document published prior to the international filing date but 
later than the priority date claimed 



later document published after the International filing date 
or priority date and not In conflict with the application but 
cited to understand the principle or theory underlying the 
Invention 

document of particular relevance; the claimed invention 
cannot be considered novel or cannot be considered to 
involve an inventive step when the document is taken alone 

document of particular relevance; the claimed invention 
cannot be considered to involve an Inventive step when the 
document is combined with one or more other such docu- 
ments, such combination being obvious to a person skilled 
in the art. 

document member of the same patent family 



Date of the actual completion of the intemational search 



30 October 2001 



Date of mailing of the international search report 



07/11/2001 



Name and mailing address of the ISA 

European Patent Office, P.B. 5818 Patentlaan 2 
NL - 2280 HV Rijswijk 
Tel. (+31-70) 340-2040. Tx. 31 651 epo nl. 
Fax: (+31-70) 340-3016 



Authorized officer 



Swart jes, H 



Form PCT/ISA/210 (second sheet) (July 1992) 



page 1 of 2 



INTERNATIONAL SEARCH REPORT ..... 

International Application No 

Pd/GB 01/03042 


C.(Contlnuation) DOCUMENTS CONSlD^^|P> BE RELEVANT ^^P' 


Category ° 


Citation of document, with indicatlon.wtiere appropriate, of tfie relevant passages ' 


Relevant to claim No. 


A 


SETTARI ANTONIN ET AL: "Advances in 
coupled geomechanical and reservoir 
modeling with applications to reservoir 
compaction" 

PROCEEDINGS OF THE 1999 15TH SYMPOSIUM ON 
RESERVOIR SIMULATION; HOUSTON, TX, USA FEB 
14-17 1999, 

1999, pages 345-357, XP001036510 
Richardson, TX, USA 

page 1, right-hand column, line 19 - line 
41 


1.7 



Form PCT/ISA/210 (continuation of second sheet) (July 1992) 



page 2 of 2 



INTERNATIONAL SEARCH REPORT 

information on patent family members 



Patent document 
cited In search report 



s 



Publication 
date 



International Application No 

PCi^GB 01/03042 



•Patent family 
merhber(s) < 



Publication 
date 



US 5417103 



23-05-1995 NONE 



Form PCT/ISA/210 (patent family amex) (July 1992) 



PATENT COOPERATION TRJKY 



From the INTERNATIONAL SEARCHING AUTHORITY 



PCT 



To: 

GECO-PRAKLA (UK) Limited 
Attn. Stool e, Brian D. 
Schlumberger House 
Buckingham Gate 
Gatwick, West Sussex RH6 ONZ 
UNITED KINGDOM 


NOTIFICATION OF TRANSMITTAL OF 
THE INTERNATIONAL SEARCH REPORT 
OR THE DECLARATION 

(PCT Rule 44.1) 


Date of mailing 

(day/month/year) 07/1 1 /2001 


Applicant's or agent's file reterence 
94.0036 


FOR FURTHER ACTION See paragraphs 1 and 4 below 


International application No. 

PCT/GB 01/03042 


International filing date 
(day/monlh/year) 06/07/2001 


Applicant 

SCHLUMBERGER EVALUATION & PRODUCTION SERVICES (UK) 



1 . [2 The applicant is hereby notified that the International Search Report has been established and is transmitted herewith. 
Filing of amendments and statement under Article 19: 

The applicant is entitled, if he so wishes, to amend the claims of the International Application (see Rule 46): 

When? The time limit for filing such amendments is normally 2 months from the date of transmittal of the 
International Search Report; however, for more details, see the notes on the accompanying sheet. 



Where? Directly to the 



International Bureau of WlPO 
34, chemin des Colomtjettes 
121 1 Geneva 20, Switzerland 
Fascimlle No.: (41-22) 740.14.35 



For more detailed instructions, see the notes on the accompanying sheet. 

2. I I The applicant is hereby notified that no International Search Report will be established and that the declaration under 
' — ' Article I7(2)(a) to that effect is transmitted herewith. 

3. With regard to the protest against payment of (an) additional fee(s) under Rule 40.2. the applicant is notified that: 

□ the protest together with the decision thereon has been transmitted to the International Bureau together with the 
applicant's request to fooA/ard the texts of both the protest and the decision thereon to the designated Offices. 

I I no decision has been made yet on the protest; the applicant will be notified as soon as a decision is made. 

4. Further action(s): The applicant is reminded of the following: 

Shortly after 18 months from the priority date, the international application will be published by the International Bureau. 
If the applicant wishes to avoid or postpone publication, a notice of withdrawal of the international application, or of the 
priority claim, must reach the International Bureau as provided in Rules QObisA and 906/S.3. respectively, before the 
completion of the technical preparations for international publication. 

Within 19 months from the priority date, a demand for international preliminary examination must be filed if the applicant 
wishes to postpone the entry into the national phase until 30 months from the priority date (in some Offices even later). 

Within 20 months from the priority date, the applicant must perform the prescribed acts for entry into the national phase 
t>efore all designated Offices which have not been elected in the demand or in a later election within 19 months from the 
priority date or could not be elected because they are not bound by Chapter 11. 
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Fax: (+31-70) 340-3016 



Authorized officer 

Jacobus Constant 



;S TO FORM PCTyiSAAi20 



These rstotes are intended to give the basto instnjdions conoeming the fiOng of amendments under aitide 19. The 
Notes are based on the requirements of the Patent Cooperation Treaty, the Regulations and the Adminbtrative Instniciions 
under that Treaty. In case of discrepancy between these Notes and those requiremento, the latter are applicable. For more 

detailed information, see also the PCT Applicant's Guide, a publication of WIPO. 

In these Notes, 'Aitide", 'Rule', and •Section' refer to the provisions of the PCT, the PCT Regulations and the PCT 
Administrative Instructions respectively. 



INSTRUCTIONS CONCERNING AMENDMENTS UNDER ARTICLE 19 

The applicant has, after having received the international search report, one opportunity to amend the claims of the 
international application. It should however be emphasized that, since all parts of the international application (claims, 
description and drawings) may be amended during the international preliminary examination procedure, there is usually 
no need to file amendments of the claims under Article 19 except where, e.g. the applicant wants the latter to be pubfahed 
for the purposes of provisional protection or has another reason for amending the claims before international pbulication. 
Furthemiore, it should be emphasized that provisional protection is available in some States only. 



What parts of ttie International application may t>e amended? 

Under Article 1 9, only the claims may be amended. 

During the international phase, the daims may also be amended (or further amended) under Article 34 before 
the International Preliminary Examining Authority. The description and drawings may only be amended under 
Artide 34 k)effore the tntemational Examining Authority. 

Upon entry into the national phase, all parts of the international application may be amended under Article 28 
or, where applicak)le, Artide 41 . 



When? Within 2 months from the date of transmittal of the tntemational search report or 16 months from the priority 

date, whichever time limit expires later. It should be noted, however, that the amendments will be considered 
as having been received on time if they are received by Ihe Intemational Bureau after the expiration of the 
applicable time limit but before the completion of the technical preparations lor intemational publication 
(Rule 46.1). 



Where not to file the amendments? 

The amendments may only be filed with the Intemational Bureau and not with the receiving Office or the 
Intemational Searching Authority (Rule 46.2). 

Where a demand for intemational preliminary examination has been^s filed, see below. 



How? Either by cancelling one or more entire daims, by adding one or more new daims or by amending the text of 

one or more of the daims as filed. 

A replacement sheet must be submitted for each sheet of the daims which, on account of an amendment or 
amendments, differa frem the sheet originally filed. 

All the daims appearing on a replacement sheet must be numbered in Arabic numerals. Where a daim is 
cancelled, no renumbering of the other claims is required. In all cases where daims are renumbered, they must 
be renumt>ered consecutively (Administrative Instructions, Section 205(b)). 

The amendments must be made In the Unguage In which thelntefnatlonal application Is to be published. 



What documents must/may accompany the amandmenU? 
Letter (Section 205(b)): 

The amendments must be sut>mttted with a letter. 

The letter will not be published with the intemational application and the amended daims. It should not be 
confused with the 'Statement under Article 1 9(1 )' (see below, under 'Statement under Artide 19(1)*). 

The letter must be In English or French, at the choice of the applicant However. If the language of the 
Intemational application Is English, the letter must be In English; If the language of the International a| 
Is French, the letter must tie In French. 
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The letter must indicate the differences between the claims as filed and the claims as amended. It must, in 
particular, indkate, in connection with each daim appearing in the international application (it being understood 
that identical indications concerning several claims may be grouped),whether 

(i) the claim is unchanged; 

(H) the claim is cancelled; 

(iii) the claim is new; 

(iv) the claim replaoes one or more claims as filed; 

(v) the claim is the result of the divbion of a claim as filed. 



The following examples illustrate the manner In which amendments must be explained in the 

accompanying letter: 

1 . [Where ofiginaUy there were 48 cl aims and after amendment of some claims there are 51 ]: 
'Oaims 1 to 29, 31 , 32, 34, 35, 37 to 48 replaced by amended claims bearing the same numbers; 
daims 30, 33 and 36 unchanged; new claims 49 to 51 added.' 

2. [Where originally there were 1 5 claims and after amendment of all claims there are 1 1 ]: 
"Claims 1 to 1 5 replaced by amended daims 1 to 1 1 

3. [Where originally there were 1 4 daims and the amendments consist in cancelling some daims and in adding 

new claims]: 

'Claims 1 to 6 and 14 unchanged; claims 7 to 13 cancelled; new claims 1 5, 16 and 17 added.* or 
*aaims 7to 13 cancelled; new claims 15, 16 and 17 added; aU other claims unchanged' 

4. pMiere various kirxis of amendments are made]: 

'Claims 1-10 unchanged; claims 1 1 to 1 3, 1 8 and 1 9 cancelled; claims 1 4, 1 5 and 1 6 replaeed by amended 
daim 14; claim 17 subdivided into an>ended daims 15, 16 and 17; new daims 20 and 21 added.' 



"Statement under article 19(1)** (Rule 46.4) 

The amendments may be accompanied by a statement explaining the amendments and indcating any impact 
that such amendments might have on the description and the drawings (which cannot be amended under 
Artide19(1)). 

The statement will be published with the international application and the amended daims. 
It must be In the language In which the International apppllcatlon la to be puMlahad. 

H must be brief, not exceeding 500 words if in English or if translated into English. 

It should not be confused with and does not replace the letter indicating the dHferenoes between the claims 
as filed and as amended, tt must be filed on a separate sheet and must t>e Identified as such by a heading, 
preferably by using the words 'Statement under Article 19(1).* 

It may not contain any disparaging comments on the international search report or the relevar>ce of citations 
contained in that report. Reference to citations, relevant to a given claim, contained in the intemational search 
report may t>e made only in connection with an amendment of that claim. 



Consequence If a demand for Intemational preliminary examination has already l>een filed 

If, at the time of filing any amendments under Article 1 9, a demand for intemational preliminary examination 
has already t>een submitted, the applicant must preferat)ly, at the same time of filing the amendments with the 
Intemational Bureau, also file a copy of such amendments with the Intemational Preliminary Examining 
Authority (see Rule 62.2(a), first sentence). 



Consequence with regard to translation of the Intemational application for entry Into the national phase 

The applicant's attention is drawn to the fact that, where upon entry into the national phase, a translation of the 
daims as amended under Artide 19 may have to be furnished to the designated/elected Offices, instead of, or 
in addition to, the translation of the daims as filed. 

For further details on the requirements of each designated/eleded Office, see Volume It of the PCT Applicant's 
Guide. 
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1 . Basis of the report 

a. With regard to the language, the international search was carried out on the basis of the international application in the 
language In which it was filed, unless othenvise indicated under this item. 



the international search was carried out on the basis of a translation of the international application furnished to this 
Authority (Rule 23.1(b)). 

With regard to any nucleotide and/or amino acid sequence disclosed in the international application, the international search 

was carried out on the basis of the sequence listing : 

I I contained in the international application in written form. 

filed together with the international application In computer readable form, 
furnished subsequently to this Authority in written form, 
furnished subsequently to this Authority In computer readble form. 

the statement that the subsequenUy furnished written sequence listing does not go beyond the disclosure in the 
international application as filed has been furnished. 

the statement that the information recorded in computer readable form is identical to the written sequence listing has been 
furnished 



□ 



□ 
□ 
□ 
□ 

□ 



2. 
3. 



I I Certain claims were found unsearchable (See Box 1). 
I I Unity of invention is laclcing (see Box II). 



4. With regard to the title, 

[X| the text is approved as submitted by the applicant. 

I I the text has been established by this Authority to read as follows: 



With regard to the abstract, 

pC] the text is approved as submitted by the applicant 



□ the text has been established, according to Rule 38.2(b). by this Authority as it appears in Box III. The applicant may. 
within one month from the date of mailing of this international search report, submit comments to this Authonty. 



The figure of the drawings to be published with the abstract is Figure No. 
I I as suggested by the applicant. 
[X| because the applicant failed to suggest a figure. 
I I because this figure better characterizes the Invention. 



li 



I I None of the figures. 
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METHOD OF DETERMINING SUBSIDENCE IN A RESERVOIR 
BACKGROUND OF THE INVENTION 

The subject matter of the present invention relates to a reservoir simulator including a 
parameter determination software which determines displacement parameters 
representing subsidence in an oilfield reservoir. 

There are many recent reports of geomechanical modelling being used predictively for 
evaluation of altemative reservoir development plans. Li the South Belridge field, 
Kem County, California, Hansen et al^ calibrated finite-element models of depletion- 
induced reservoir compaction and surface subsidence with observed measurements. 
The stress model was then used predictively to develop strategies to minimise 
additional subsidence and fissuring as well as to reduce axial compressive type casing 
damage. Berumen et afi developed an overall geomechanical model of the Wilcox 
sands in the Arcabuz-Culebra field in the Burgos Basin, northern Mexico. This 
model, combined with hydraulic firacture mapping together with fracture and reservoir 
engineering studies, was used to optimise firacture treatment designs and improve the 
plarming of well location and spacing. 

The subject of fluid flow equations which are solved together with rock force balance 
equations has been discussed extensively in the literature. Kojic and Cheatham^'"^ 
present a lucid treatment of the theory of plasticity of porous media with fluid flow. 
Both the elastic and plastic deformation of the porous medium containing a moving 
fluid is analyzed as a motion of a solid-fluid mixture. Corapcioglu and Bear^ present 
an early review of land subsidence modelling and then present a model of land 
subsidence as a result of pumping firom an artesian aquifer. 
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Demirdzic et af'^ have advocated the use of finite volume methods for numerical 
solution of the stress equations both in complex domains as well as for thermo-elastic- 
plastic problems. 

A coupling of a conventional stress-analysis code with a standard finite difference 
reservoir simulator is outlined by Settari and Walters^. The term "partial coupling" is 
used because the rock stress and flow equations are solved separately for each time 
increment. Pressure and temperature changes as calculated in the reservoir simulator 
are passed to the geomechanical simulator. Updated strains and stresses are passed to 
the reservoir simulator which then computes porosity and permeability. Issues such as 
sand production, subsidence, compaction that influence rock mass conservation are 
handled in the stress-analysis code. This method will solve the problem as rigorously 
as a fully coupled (simultaneous) solution if iterated to fiill convergence. An explicit 
coupling, i.e. a single iteration of the stress model, is advocated for computational 
efficiency. 

The use of a finite element stress simulator with a coupled fluid flow option is 
discussed by Heffer et af and by Gutierrez and Lewis^^. 

Standard conraiercial reservoir simulators use a single scalar parameter, the pore 
compressibility, as discussed by Geertsma^^ to account for the pressure changes due to 
volumetric changes in the rock. These codes generally allow permeabiUty to be 
modified as a function of pore pressure through a table. This approach will not be 
adequate when the flow parameters exhibit a significant variation with rock stress. 
Holt^^ found that for a weak sandstone, permeability reduction was more pronoimced 
under non-hydrostatic applied stress, compared with the slight decrease measured 
under hydrostatic loading. Rhett and Teufel^^ have shown a rapid decline in matrix 
permeability with increase in effective stress. Ferfera et al^^ worked with a 20% 
porosity sandstone and found permeability reductions as high as 60% depending on 
the relative influence of the mean effective stress and the differential stress. Teufel et 
al^^ and Teufel and Rhett^^ foimd, contrary to the assumption that permeabiUty will 
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decrease with reservoir compaction and porosity reduction, that shear failure had a 
beneficial influence on production through an increase in the fracture density. 

SUMMARY OF THE INVENTION 

Oil recovery operations are seeing increased use of integrated geomechanical and 
reservoir engineering to help manage fields. This trend is partly a result of newer, 
more sophisticated measurements that are demonstrating that variations in reservoir 
deliverability are related to interactions between changing fluid pressures, rock 
stresses and flow parameters such as permeability. Several recent studies, for 
example, have used finite-element models of the rock stress to complement the 
standard reservoir simulation. 

This specification discusses current work pertaining to: fully and partially coupling 
geomechanical elastic/plastic rock stress equations to a conmiercial reservoir 
simulator. This finite difference simulator has black-oil, compositional and thermal 
modes and all of these are available with the geomechanics option. In this 
specification, the implementation of the aforementioned stress equations into the code, 
hereinafter called the 'Parameter Determination Software*, is discussed. Some work 
on benchmarking against an industry standard stress code is also shown as well as an 
example of the coupled stress/fluid flow. The goal in developing this technology 
within the simvdator is to provide a stable, comprehensive geomechanical option that 
is practical for large-scale reservoir simulation. 

Further scope of applicability of the present invention will become apparent from the 
detailed description presented hereinafter. It should be understood, however, that the 
detailed description and the specific examples, while representing a preferred 
embodiment of the present invention, are given by way of illustration only, since 
various changes and modifications within the spirit and scope of the invention will 
become obvious to one skilled in the art from a reading of the following detailed 
description. 
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BRIEF DESCRIPTION OF THE DRAWINGS 

A full understanding of the present invention will be obtained from the detailed 
description of the preferred embodiment presented hereinbelow, and the 
accompanying drawings, which are given by way of illustration only and are not 
intended to be limitative of the present invention, and wherein: 

figure 1 illustrates control volume (light grey grid) for Rock Force Balance; 

figure 2 illustrates comparison of Rock Displacement Predictions (case 1); 

figures 3A and 3B illustrate comparison of Rock Displacement Predictions (case 2); 

figures 4al, 4a2, 4b 1, 4b2, 4c 1, and 4c2 illustrate total normal stress in Mid X-Y 
plane; 

figure 4d illustrates Table 1, and Table 1 further illustrates simulation parameters used 
in comparing ABAQUS and Reservoir Simulator Predictions); 

figures 5-9 illustrate examples of stractured and unstructured grids; 

figures 10 and 11a illustrate examples of a staggered grid; 

figures lib, 12 and 13 illustrate the 'Tlogrid" software and the "EcUpse simulator 
software" and the generation of 'simulation results' which are displayed on a 3-D 
viewer; 

figure 14 illustrates an example of the 'simulation results'; 

figure 15 illustrates the fact that the "Eclipse simulator software" includes the 
'Parameter Determination Software' of the present invention; 
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figure 16 illustrates the results which are generated when the 'Parameter 
Determination Software' is executed by a processor of a workstation, the results being 
presented and displayed on a recorder or display device; 

figure 17 illustrates how the results generated in figure 16 assist in determining 
subsidence in an oilfield reservoir; and 

figure 18 is a detailed flowchart of the Tarameter Determination Software' of the 
present invention. 



DESCRIPTION OF THE INVENTION 



Elastic Stress Equations 

Steady state "rock momentum balance equations" in the x, y and z directions can be 
written as^^: 



dx dy 

3(T„ 3r_ 3t_, ^ „ ^ 
— —+F +P =0. 

dz ' " 



dy dx 



.(1) 



dz dx 



dy 



Here F is the body force, P is the interaction force between the solid and the fluid. 
This interaction force is discussed further below. 
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The elastic normal stresses o and shear stresses t can be expressed in terms of strains, 
£ and y , as: 



Constants G, also known as the modulus of rigidity, and A are Lamp's constants. 
They are functions of Young's modulus, E, and Poisson's ratio, v . Strains e^^ y^are 

defined in terms of displacements in the x,y,z directions, namely u,v and w. Thus 



e —^^ 

dy 
_ dw 

'''^ (3) 
_ dv dw 



_ dw du 



Within the simulator, the rock force is calculated in units of Newtons, Ibf, or dynes 
depending on whether the user has chosen mks, field or cgs units for the simvilation. 
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Gridding 

Before returning to a discussion of the fluid component mass conservation and rock 
mass conservation equations, the special gridding for the rock force (momentum) 
balance equations is briefly outlined. 

To implement these rock momentum balance equations in a finite difference 
simulator, separate control volumes for balancing the force are used in each 
orthogonal grid direction. These control volumes are staggered in each of the 
coordinate directions. 

In Figure 1, an x-direction rock force balance control volume is shown. Also shown 
in this Figure is the finite difference discretization used. The grid with the darker 
black lines in this figure is the standard reservoir simulation grid which are control 
volumes for fluid component and rock mass conservation. Rock displacements u, v 
and w are defined on the faces of this grid. Rock momentum control volumes in a 
particular coordinate direction are centred on the rock displacement in that direction. 
In Figure 1 is shown the x-momentum control volume (light grey grid) in the centre of 
which is the x-direction rock displacement, u, which is defined on the face between 
two control volxmies for mass conservation (black grid). Similar rock momentum 
control volumes are set up in the y and z directions around v and w, the y and z rock 
displacements 

This grid is fixed in time. Rock will flow elastically or plastically in the mass 
conservation grid. A correction term to Darcy's law that accounts for the extra 
transport of fluid component mass into an adjoining mass conservation ceU due to 
rock movement into that cell can be calculated in the simulator and is discussed 
below. 

We have found the staggered gridding to be the most natural for differencing purposes 
and for setting of boundary conditions. It is also the most accurate because it achieves 
second order accuracy when calculating the normal stresses and mixed first/second 
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order accuracy when calculating the shear stresses. To illustrate, if the user has input 
the noLore general stiffness matrix^^, 5^^, instead of the Young's moduli and Poisson 

ratios, then the differenced normal stress on the x+ side of the rock momentum 
balance control volume in Figure 1 can be written 



Ax. . , 



2x+ 



k 



Another consequence of the staggered gridding is that the force balance equations (1) 
give rise to a tfairteen-banded jacobian for three-dimensional calculations, hi addition 
to the standard diagonal and six off-diagonal bands, there are six more connecting 
points in the differencing stencil which contribute the other bands. 

Fluid Component Conservation, Rock Mass Conservation and Volume Balance 
Equations 

Compositional and thermal simulation require the user to specify the number of fluid 
components. Fluid components are conserved in the standard user-specified grid 
blocks as noted above. Rock also flows either elastically or plastically, in addition to 
fluid, amongst these grid cells. A rock mass conservation equation can be written in 
pseudo-differenced form as: 

v/>„.i(i-<.r-(i-</>)'j= 

faces 



where 
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V = bulk volume of the rock mass conservation cell 
Prodfc = specific density of rock 
^ = porosity of rock (variable) 
u,v,w = x,y,z total rock displacement 

Superscripts n+1, n refer to the time level - n refers to the last time step, n+1 refers to 
the current or latest time level. The superscript "up" refers to the upstream of the rock 
displacement change over the time step. Flow of rock into an adjoining cell originates 
from the cell away from which rock flows during the timestep, as denoted by 
M,v,w^^^ - w,v,w^^^ in the above equation (5). Total rock displacement is the sum of 

elastic and plastic rock displacement. The elastic rock displacements in the 
coordinate directions and porosity are included in the simulator variable set. Rock 
properties, for example relative permeability endpoints, are convected. 

Fluid component mass conservation assumes the usual form and is discussed 
elsewhere^^*^^. 

In order that the fluid volumes fit into the available pore volume, a final fluid phase 
volume balance equation is included. It is often written as the sum of fluid phase 
saturations (volume of fluid divided by pore volume) equal to unity. 

Additional Enhancements to the Balance Equations 

A term that describes the change of rock momentum in time, namely vi\^^ , where 

wz^^is the mass of rock in the rock balance grid cell, can also optionally be switched 

on by the user. By default it is onadtted because this term is small and makes litde 
contribution to the force balance. Also, when it is not specified, wave solutions of the 
form f(x±ct) are not calculated where c is some characteristic rock 
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compression/shear velocity. These wavelets may cause instability and usually are of 
little interest over standard simulation time scales. 

Secondly, as discussed by Corapcioglu and Bear^, in a deforming porous medium it is 
the specific discharge relative to the moving solid, q^, that is expressed by Darcy's 
law 

qr=q-<pV,^-KVH (6) 

where 

q = specific discharge of a fluid component (flux) 
^ = porosity 

= velocity of the solid 

K = fluid phase permeability 

H = pressure head, including the gravitational head 

This term is, again, omitted by default but can be included if desired. In problems 
where the elastic and plastic rock flow is small and near steady state, there is not much 
contribution from this term. In situations where the pressure field is changing 
continuously and strongly, for example in a thermal cyclic huff and puff scenario 
where the near wellbore pressure can vary by more than 500 psi through a short cycle 
time period, this term can have a strong effect on the simulation. In particular, it can 
reduce the unphysically high injection pressures normally predicted by a conventional 
simulator. There is an additional CPU expense associated with this term since it can 
add one or more newton iterations to a timestep. 

Interaction between the fluid and the rock 

The interactive force between the fluid and solid, denoted Pin equation (1), can be 
expressed in the form^ 

10 



SUBSTITUTE SHEET (RULE 26) 



wo 02/06857 



PCT/GBOl/03042 



(7) 



where 



inertial force of the fluid per unit bulk volume 



body force acting on the fluid per unit fluid volume 



stress tensor of the fluid 



The stress tensor of the fluid can be related to the pore fluid pressure, following 
Terzaghi's^^ concept of intergranular stress, 



where / is the unit tensor. We neglect the effect of the inertial forces of the fluid and 
write the interaction force as 



The pore pressures are calculated at the centres of the mass conservation cells so that 
the pore pressure gradient in equation (9), Vp , is centred in the rock momentum 
conservation control volume as it should be. 

Interaction of the rock stress on flow parameters, in particular fluid permeability, is a 
subject that is the topic of much current research. As noted above in the introduction, 
many authors are now carrying out triaxial compression tests on core and rock 
samples in which fluid is flowing in the stressed rock matrix. Most conunercial 
reservoir simulators contain permeability or transmissibility modifiers as a function 
of fluid pore pressure. This feature was used in a geomechanics simulation study by 




pi. 



(8) 



(9) 
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Davies and Davies where the relationship between stress and permeability was 
calculated based on the relationship between the stress and porosity, and between 
permeability and porosity for various rock types. Yale and Crawford"^ discuss the 
work of Holt^^, Teufel and Rhett^^*^^'^^ and others. They also used critical state theory 
to model observed stress-permeability experimental data. It was concluded that this 
theory did well represent the observed data. 

We are currently exploring several constitutive relationships between stress and 
permeability. A mean effective stress versus permeability is one of these 
relationships, another is the classical permeability-porosity. Another possibility is 
permeability versus 0,-03, the difference between the maximum and nndnimum 

principal stresses, as displayed on the web by GFZ-Potsdam^"^. The success of the 
critical state theory as discussed in the above paragraph is also of interest. 

Variables and Solution of the Equations 

We solve for moles of all fluid components per unit pore volume, energy per unit bulk 
volume if thermal simulation is specified, pressure of one of the fluid phases, elastic 
rock displacements in the x,y,z directions and porosity. Some of these variables are 
eliminated before entering the solver. In particular, one of the fluid molar densities is 
eliminated using the volume balance equation. Usually, a fully implicit solution of all 
equations is unnecessary because in many parts of the grid the fluid throughput in a 
mass conservation grid cell is not high. In this case an adaptive implicit solution of 
the equation set is used in which some of the grid blocks are solved fiilly implicitiy, 
while others are solved using the IMPES (implicit pressure, explicit saturation) 
solution scheme. Often the IMPES scheme is used everywhere although the 
timestepping of this particular scheme is then limited by the throughput. The rock 
force balance equations are always solved fully implicitly, even in grid cells which are 
IMPES. 

The linear solver used to solve the coupled system of non-linear equations is a nested 
factorization technique^. One such linear solver is disclosed in US Patent 6,230,101 
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to Wallis, the disclosure of which is incorporated by reference into this specification. 
Although the rock stress equations (1) are nonlinear because of the impUcit porosity 
and highly elliptic and, in addition, are placed next to mass conservation equations 
which display both hyperbohc and parabolic form in the jacobian, the solver has 
demonstrated good robustness in solving the coupled system of equations. 

Although some authors^'^^ have described fiilly or partially coupled schemes where 
porosity is calculated directly using rock dilation, f^H-fy + f^, we have chosen the 

additional rock conservation equation for two reasons. Firstly, the approach is more 
general and allows rock to be produced in cases of sand fluidization, wellbore 
instability, etc. Secondly, as discussed further below, at the beginning of a new 
timestep a large plastic displacement of the rock is accounted for exactly by including 
it in the implicit rock mass balance for this timestep. 

Plasticity 

An explicit plastic calculation has been implemented. The Mohr-Coulomb plastic 
calculation is available. The user can specify the cohesion and angle of internal 
friction in chosen regions in the grid. At the start of a new timestep, principal stresses 
are calculated in each rock force-balance control volume based on the stress field from 
the last timestep. Diameters of the smallest and largest Mohr circles are calculated 
and the test is made as to whether failure has occurred. If so, the displacement around 
which the control volume is centred is altered to bring the circle back to the failure 
envelope. This corrective displacement is stored as a plastic displacement. These 
plastic displacements, which can be referred to as (Up, Vp, Wp), are part of the total 
displacements (u, v, w) which in turn are a sum of the elastic and plastic contributions. 
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Boundary Conditions and Initialisation 

Either stress or displacement boundary conditions can be chosen. Each side of the 
simulation grid, that is x-, x+, y-, y+ and z-, z+, can have one of the above boundary 
conditions prescribed. Usually the bottom of the grid will have a zero displacement 
condition while the sides and top will have specified stresses. 

Initialisation of a reservoir with stress boundary conditions requires an initial small 
simulation to equilibrate the simulated field. Before the wells are turned on, boundary 
stresses are ramped up at intervals until the desired vertical and horizontal stresses are 
achieved. The pore pressure is maintained at an initial level by including an extra 
production well with a BHP limit set to this level. Compression of the reservoir by 
the boundary stresses causes the rock matrix to compress which in tum forces some of 
the fluid out of this well. Huid-in-place must be checked before starting the 
simulation. 

User Interface 

The user enters a Young's modulus, Poisson ratio and specific rock density for each 
grid block. It is also possible to input the full stiffness matrix, as illustrated above in 
equation (4), to allow an anisotropic stress tensor. If a plastic calculation is required, 
the user enters a cohesion and angle of internal friction in different regions of the grid. 

Efficiency of Geomechanical Computations 

Because there are several extra equations to be set up and solved, there is a sizable 
overhead associated with the geomechanics calculation. This overhead is offset by 
several factors. Firsfly, the coupled set of equations is stable and robust. Secondly, 
implementation of the jacobian and right-hand-side setup includes vectorization 
wherever possible. Also, parallel options exist within the simulator. 
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On average, for three dimensional IMPES sinniulations we have noted that run times 
are about three times as long as the run time needed for the same calculation without 
geomechanics. This is efficient considering the IMPES calculation only solves a 
single equation, whereas four additional equations are solved when geomechanics is 
included. Both the efficiency of setting up the equations and the robustness of the 
solver contribute to the overall efficiency. The simulator is capable of achieving good 
parallel speedups and this also helps to offset the increased computational demands. 

Benchmarking Against an Industry Standard Stress Code 

We have chosen the ABAQUS^^ stress simulator to test the rock stress calculations in 
our reservoir simulator. ABAQUS was chosen for two reasons. Firstly, it is widely 
accepted and established. Secondly, it contains a simple, single-phase flow-in-porous- 
media option. 

Two elastic, one-dimensional problenos were designed and run in both simulators. The 
test cases differed only in that one featured an average sandstone, the other a weak 
sandstone. They were designed to have large pore pressure gradients so that the fluid- 
rock interaction was pronounced. 

In figure 4d, Table 1 presents the simulation parameters used in the test problems. 
These are listed in the table as Cases 1 and 2. A water injection rate of 1000 BBL/D 
was used. The injector was situated in one comer, the producer in the far comer. 
Boundaries were rigid. Dimensions of the grid were chosen as representative of a 
small pattern of injectors and producers. The reservoir simulator was ran to steady 
state for these cases. In figure 4d, under Table 1, in all cases, porosity = 0.33, fluid 
permeability = 6mD, initial pressure = 4790 psia, injection/production rate = 1000 
BBL/D and grid dimensions were 1000 ft x 1000 ft x 20f. 

In Figure 2, figure 2 shows a comparison of rock displacements as predicted by the 
two simulators for Case 1, the average strength sandstone. 
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In Figures 3 A and 3B, these figures illustrate a comparison for the weak rock. Much 
larger rock displacements are evident as a result of the very low Young's modulus. In 
both cases, the comparison against AB AQUS was within engineering accuracy. 

Other Test Cases 

In figure 4d, Table 1 is illustrated (Table 1 illustrating simulation parameters used in 
comparing ABAQUS and Reservoir Simulator Predictions), We have also run 
ABAQUS on a further set of 4 problems. Table 1, shown in figure 4d, lists these as 
Cases 3 to 6. These are a pair of two-dimensional and a pair of three-dimensional 
runs. For equilibration purposes, the ABAQUS simulations used boundary conditions 
that are not available in our reservoir simulator. We are currently working to remedy 
this. Initial comparisons in which we have approximated the ABAQUS simulation 
have been encouraging. Both quantitative and qualitative agreement has been 
achieved. 

A plastic case was run in ABAQUS to compare against the reservoir simulator plastic 
calculation. It is the same as Case 6 in Table 1 of figure 4d. A cohesion and angle of 
intemal friction was input to define the Mohr-Coulomb failure curve. Differences in 
boundary conditions in the two simulators, as discussed above, have prohibited a 
direct comparison but, again, there is qualitative agreement when comparing plastic 
displacements. The plastic calculation is robust and timestepping is not adversely 
affected. 

In Figures 4al and 4a2, 4b 1 and 4b2, and 4c 1 and 4c2, we have simulated some larger 
fields with the geomechanics option. An example is shown in figures 4al and 4a2. 
In this three-dimensional, compositional (10 component) problem, 7 injection and 6 
production wells contribute to the rock stress. Contours of the total normal stress, 
o^'irOy'^o^ which is also the first stress invariant, are shown in the mid-xy plane of 

the simulation. This first stress invariant is important because it largely governs 
compaction processes. Timestepping with the geomechanics calculation included was 
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the same as without There was a factor of 2.8 slowdown in the overall runtime with 
geomechanics. 

In summary, the implementation of elastic/plastic stress equations into a commercial, 
finite difference reservoir simulator has been discussed in this specification. Our goal 
of demonstrating robustness together with a comprehensive geomechanical option has 
been achieved. Work is continuing to develop this option further. In particular: 

1. The coupled system of equations has demonstrated good stability and robustness. 
In some test cases, there was no difference in timestepping when the simulation 
was run with or without the geomechanics option. 

2. Our current capabilities are comprehensive enough to predict both elastic and 
plastic rock displacement. The basic design is general enough to allow future 
engineering development. This could include wellbore stability and failure, sand 
production, more accurate fault modelling. 

3. We have indicated at various points of this paper our areas of current work and 
investigation. In particular these include stress-permeability relationships and 
modifying our geomechanical boundary conditions to compare more exactly with 
ABAQUS. We are also currendy working with a large-scale simulation of a very 
brittle carbonate field to predict subsidence. 

Nomenclature 



E 


= Young's modulus 


G 


= modulus of rigidity. Lame constant 


H 


= pressure head, including the gravitational head 


K 


= fluid phase permeability of rock matrix 


q 


= single phase fluid flux in a porous medium 


u 


= rock displacement in the x-direction 


V 


= rock displacement in the y-direction 




= velocity of the solid, as related to Darcy's law 


W 


= rock displacement in the z-direction 
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o^ y ^ = elastic nonnal rock stress in x,y,z directions 
= x,y,2 elongation strains 

Jxy.yz,zx = shcaT Strains 

A = Lame constant 

(p = porosity 

V = Poisson's ratio 

Pnck = ^^^^ specific density 

^xy,yz,xz ~ ^l^^c shcaT strcss 

DETAILED DESCRIPTION OF THE INVENTION 

Referring to figures 5 and 6, structured and unstructured grids are illustrated. 

In figure S, an earth formation IS is illustrated, the formation 15 including four (4) 
horizons 13 which traverse the longitudinal extent of the formation 15 in figure 5. 
Recall that a "horizon" 13 is defined to be the top surface of an earth formation layer, 
the earth formation layer comprising, for example, sand or shale or limestone, etc. 

A "grid" is located intermediate to the horizon layers 13. More particularly, a "grid" 
is formed: (1) in between the horizons 13, (2) on top of the uppermost horizon 13, and 
(3) below the lowermost horizon 13 in the Earth formation 15. When gridding the 
formation 15, the act or function of "gridding" the formation 15 includes the step of 
dividing the formation 15 into a multitude of individual "cells" which, when 
coimected together, comprise the "grid". 

In figure 6, for example, the Earth formation 15 includes an uppermost horizon 13a 
and a lowermost horizon 13b which is separated from the uppermost horizon 13a by 
an intermediate earth formation layer 15a, The intermediate earth formation layer 15a 
includes, for example, a sand layer or a shale layer or a limestone layer, etc. A 
particular "gridding software" will "grid" the earth formation layer 15a; that is, the 
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formation layer 15a will be divided up, by the "gridding software" into a multitude of 
cells 15al. 

In the prior art, a "gridding software*' product known as "Grid" was marketed by 
GeoQuest, a division of Schlumberger Technology Corporation, Abingdon, the United 
Kingdom (UK). The "Grid" software would divide the formation layers 15a into a 
multitude of cells. However, each of the multitude of cells were approximately 
"rectangular" in cross sectional shape. 

In addition, another "gridding software", known as "Petragrid" is disclosed in U.S. 
Patent 6,018,497, the disclosure of which is incorporated by reference into this 
specification. The "Petragrid" gridding software will 'grid' the formation with 
triangularly/tetrahedrally shaped grid cells, known as "unstructured grid cells". 

In addition, another "gridding software", known as "Hogrid", is disclosed in U.S. 
Patent 6,106,561, the disclosure of which is incorporated by reference into this 
specification. The "Flogrid" gridding software includes the "Petragrid" gridding 
software; however, in addition, the "Flogrid" gridding software will *grid' the 
formation with rectangularly shaped grid cells, known as **structured grid cells". 

In figure 6, the cells 15al are shown to be approximately "rectangular" in cross 
sectional shape. These grid cells are "stractured grid cells". 

In figure 5, however, a multitude of grid cells 15al have been formed in the earth 
formation 15 intermediate the horizons 13, and each cell 15al has a cross sectional 
shape that, in addition to being approximately '^rectangular^' in cross section, is either 
approximately "polygonal" or "tetrahedral" in cross section. Figure 5 clearly shows a 
multitude of cells 15al where each cell 15al has a cross sectional shape which is 
either approximately "polygonal" or "tetrahedral" in cross sectional shape (i.e., an 
"unstructured grid cell") in addition to being approximately "rectangular" in cross 
sectional shape 
(i.e., a "structured grid ceU"). 
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Referring to figures 7, 8, and 9, other examples of "structured'* and "unstructured" 
grid cells are illustrated. Figure 7 illustrates examples of "un-structured" grid cells 
having a triangular/tetrahedral cross sectional shape. Figtire 8 illustrates examples of 
"structured" grid cells having an approximate "rectangular" cross sectional shape. 
Figure 9 illustrates further examples of "unstructured" grid cells having a 
triangular/tetrahedral cross sectional shape. 

Referring to figures 10 and 11a, a "staggered grid" will now be defined with reference 
to figures 10 and 11a. 

In figure 10, a plurality of "structured" grid cells are illustrated (i.e., cells having an 
approximately rectangularly shaped cross sectional shape). However, in figure 10, a 
pair of "staggered grid cells" are also illustrated. Assuming that a 'first stmctured or 
unstmctured grid cell' is disposed directly adjacent a 'second neighboring structured 
or vmstmctured grid cell*, a "staggered grid cell" can be defined as one which consists 
of a *first half and a 'second half, where the *first half of the 'staggered grid cell* 
comprises one-half of the 'first structured or unstructured grid cell* and the 'second 
half of the 'staggered grid cell* comprises one-half of the 'second neighboring 
stmctured or unstmctured grid cell*. 

For example, in figure 10, the staggered grid cell 20 includes the one-half 20a of a 
first structured grid cell and the one-half 20b of the neighboring, second stmctured 
grid cell. Similarly, the staggered grid cell 22 includes the one-half 22a of a furst 
stmctured grid cell and the one-half 22b of a neighboring, second stmctured grid cell. 

In figxire 1 la, a pair of "unstmctured** grid cells are illustrated (i.e., cells each having 
an approximately triangularly or tetrahedraUy shaped cross sectional shape). In figure 
1 la, a staggered grid cell 24 consists of a first half 24a from one unstmctured grid cell 
26 and a second half 24b from a second unstmctured grid cell 28. 
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Referring to figures 1 lb, 12, 13, and 14, a computer workstation 44 or other computer 
system 44 is illustrated. 

In figure lib, the computer workstation 44 in figure lib includes a system bus, a 
workstation processor 44d connected to the bus, a workstation memory 44a connected 
to the bus, and a recorder or display or 3D viewer 44e connected to the bus. A CD« 
Rom 46 is inserted into the workstation and the following software is loaded from the 
CD-Rom 46 and into the workstation memory 44a: (1) a 'Tlogrid" software 46a, and 
(2) an Eclipse simulator software 46b. The "Flogrid" software is described in U.S. 
Patent 6,106,561 to Farmer, the disclosure of which has already been incorporated by 
reference into this specification. Input data is provided to and input to the workstation 
44: (1) a well log output record 42, and (2) a reduced seismic data output record 24a. 

In figure 12, the workstation memory 44a of figure lib stores the Hogrid software 
46a and the Eclipse simulator software 46b. The Flogrid software 46a includes a 
reservoir data store, a reservoir framework, a stmctured gridder, a "Petragrid" Un- 
structured gridder, and an Upscaler. The outputs 47 from the Upscaler and from the 
Petragrid Unstractured gridder are provided as inputs to the Eclipse simulator 
software 46b. The "Petragrid" Un-structured gridder is disclosed in U.S. Patent 
6,018,497 to Gunasekera, the disclosure of which has already been incorporated by 
reference into this specification. The Eclipse simulator software 46b, responsive to 
the outputs 47, generates a set of 'simulation results' 48 which are provided to the 3D 
viewer 44e. The Eclipse simulator software 46b includes a "Linear Solver" which is 
discussed in U.S. Patent 6,230,101 to Wallis, the disclosure of which is incorporated 
by reference into this specification. 

In figure 13, to simmiarize figures lib and 12, output data 47 is generated from the 
Upscaler and the Petragrid un-structured gridder in the Flogrid software 46a, and that 
output data 47 is provided to the Eclipse simulator software 46b, along with the 
outputs of other programs 49. Responsive thereto, the Eclipse simulator software 46b 
generates a set of simulation results 48 which are, in tum, provided to the 3D viewer 
44e. 
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In figure 14, one example of the simulation results 48, which are displayed on the 3D 
viewer 44e of figure 13, is illustrated. 

Referring to figure 15, as previously discussed, the workstation memory 44a stores the 
Eclipse simulator software 46b. However, in accordance with one aspect of the 
present invention, the Eclipse simulator software 46b includes a Parameter 
Determination software 50. The "parameters" are determined for each "staggered grid 
cell", such as each of the "staggered grid cells" shown in figure 10 and 11a of the 
drawings. The "parameters", which are determined by the Parameter Determination 
software 50 for each * staggered grid cell', include the following "parameters", which 
have been previously discussed in the "Description of the Invention" portion of this 
specification: 

u = rock displacement in the x-direction 
V = rock displacement in the y-direction 
w = rock displacement in the z-direction 

Therefore, the parameter (u, v, w) represent the movement of the rock. When the 
aforementioned parameters (u, v, w) have been determined, the following "additional 
parameters" are determined for each 'staggered grid cell' from and in response to the 
parameters (u, v, w), and these "additional parameters" have also been discussed in 
the "Description of the Invention" portion of this specification: 
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^jc.y.z = x,y,z elongation Strains 
7jcy.yz,zx = shcaT stiains 

o^ y^^ = elastic normal rock Stress in x,y,z directions 

^xy,yz.xz " clastic shear stress 

^ = porosity of rock (variable) 

Up = plastic rock displacement in the x-direction 

Vp = plastic rock displacement in the y-direction 

Wp = plastic rock displacement in the z-direction 

Recall equation (3) set forth in the "Description of the Invention" portion of this 
specification, as follows: 
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(3) 



We start by estimating values for (u, v, w), which is hereinafter referred to as the 
'estimating procedure'. Therefore, when (u, v, w) is estimated by the Parameter 
Detemiination software 50 of the present invention, the values of "f^^y^" (the 'x,y,z 

elongation strains') and ry^yi^Tx"' *shear strains') are detemained from the above 

equation (3). When and "^j^^^^^^^" are determined from equation (3), the 

values of "0^,^ 'elastic normal rock stress in x,y,z directions') and "T^^y y^ j^^" 
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(the *elastic shear stress') are further determined from the above described equation 
(2), as follows: 



Now that "Oj^y^" and "t^^.^ ^^" are known, the previously described set of *rock 
momentum balance' differential equations [that is, equation (1)] are solved using the 
aforementioned known values of "cJ^.^ "'^xyo^.xz*'- After the 'rock momentum 

balance' differential equations are solved, if the resultant 'residuals' are determined to 
be approximately equal to zero (0), the previously estimated values of (u, v, w) are 
thereby determined to represent 'accurate rock displacement parameters' for the 
reservoir. Subsidence is determined from the 'accurate rock displacement parameters' 
by solving a final set of failure criterion equations, which comprise the residuals and 
any derivatives, to determine a set of rock plastic displacements (Up = plastic rock 
displacement in the x-direction, Vp = plastic rock displacement in the y-direction, and 
Wp = plastic rock displacement in the z-direction) forming a part of the rock 
displacement parameters 



The Parameter Determination software 50 will also determine the porosity of 
rock, which is a variable set forth in equation (5) set forth above. 

Referring to figure 16, when the parameters (u, v, w, e^ ^^^ , J:^.^y^,^. ^xy,yz.xz- 
^ , Up, Vp, Wp) have been determined, as discussed above, for each staggered grid block 
or cell, the Recorder or Display or 3D Viewer 44e of figure lib will illustrate a 
'Table" similar to the 'Table" shown in iBgure 16. 




zx 



(U, V, W). 
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In the Table of figure 16, for staggered grid block 1 ( "SGBl"), a first set of the 
parameters "(u, v, w, e^ ^ ^ , J^^y^^^, <^x.y,z^ ^xy,yz.:cz> ^ ' ^P' ^p* ^p)" will be generated 
and displayed on the Recorder or Display device 44e. 

Similarly, in the Table of figure 16, for staggered grid block 2 ("SGB2")5 a second set 
of the parameters "(u, v, w, £^^^, o^^,,, ^xy.yz.xz^ ^P' V Wp)" will be 

displayed on the Recorder or Display device 44e. 

Similarly, in the Table of figure 16, for staggered grid block "n" ("SGBn"), an n-th set 
of the parameters "(u, v, w, e^^^, Jxy.yz^zx^ <^x.y.z» '^xy.yz.xz^ ^» ^p» V Wp)" will be 
generated and displayed on the Recorder or Display device 44e. 

Referring to figure 17, when the "Table" of figure 16 is generated and displayed on 
the display device 44e, it is now possible to determine 'Subsidence' in an oilfield 
reservoir. The term 'subsidence' refers to a failing or 'giving away' of the seabed 
floor. In figure 17, a drilling rig 52 is situated at the surface of the sea 54. The 
drilling rig 52 withdraws underground deposits of hydrocarbon (e.g., oil) and other 
fluids (e.g. water) from an undersea reservoir 56 which is located below the seabed 
floor 58. Over a period of, for example, ten years, a certain amoimt of subsidence 60 
occurs due to the withdrawal of the underground deposits of hydrocarbon and water 
from the undersea reservoir 56, the subsidence 60 producing a lowering of the seabed 
floor 58 a certain amount which corresponds to the subsidence 60. 

In figure 17, the Parameter Determination Software 50 of figure 15 determines the 
aforementioned subsidence 60 by fibrst determining (via the 'estimating procedure' 
discussed above) the * accurate rock displacement parameters (u, v, w)' for the 
reservoir. These 'accurate rock displacement parameters (u, v, w)' represent the 
'movement of the rock', or, in our example, the displacements of the seabed floor 58 
in the (x, y, z) directions over the ten year example period. These displacements of 
the seabed floor 58 occur as a result of the subsidence 60 (that is, the 'failing* or the 
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'giving away' of the ground or seabed floor) illustrated in figure 17. Therefore, the 
parameters (u, v, w) represent or characterize 'rock movement in the x, y, and 2 
directions' which, in turn, represent or characterize the *subsidence 60'. Note that the 
'accurate rock displacement parameters (u, v, w)' may include: (1) a displacement that 
is 'elastic', and (2) a displacement that is 'plastic' which is denoted (up, Vp, Wp). 
These plastic displacements (Up, Vp, Wp) were referred to above under 'Description of 
the Invention' in the section entitied 'Plasticity'. 'Subsidence* is the result of a 
'failure' of the rock, for example by crushing, cracking or some other failure 
mechanism. When the rock has failed, some of its displacement is not recoverable 
when tihie original conditions are imposed, and it is the presence of this unrecoverable 
displacement (Up, Vp, Wp) that characterizes the 'subsidence' 60. Also note that the 
undersea formation in figure 17 is gridded with structured and unstructured grids, and, 
as a result, the undersea formation in figure 17 is gridded with "staggered grids", as 
graphically illustrated by the staggered grid 62 in figure 17. The flowchart of figure 18 
will describe the method, practiced by the Parameter Determination Software 50 of 
figure 15, for determining the 'accurate rock displacement parameters (u, v, w)' which 
represent these displacements in the (x, y, z) directions due to the subsidence 60. 

Referring to figure 18, a flowchart is illustrated which describes the method practiced 
by the Parameter Determination Software 50 for determining the parameters (u, v, w) 
representing the rock displacements in the (x, y, z) directions (i.e., the 'movement of 
the rock') due to the subsidence illustrated in figure 17. 

In figure 18, the Parameter Determination Software 50 will determine the 'rock 
movement' parameters (u, v, w) in accordance with the following steps. With 
reference to figure 18, each step of the Parameter Determination software 50 wiU be 
discussed in detail below as follows: 
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1. Input user defined rock mechanical properties - block 50a 

The user inputs a lot of data, called the "input file", comprising ke3rwords that define 
the particular simulation that is desired. Part of these keywords represent the 
mechanical properties of the rock. These "mechanical properties of the rock" are the 
kind of properties that are needed to drive the equations set forth above known as the 
"rock momentum balance equations" which are identified above as "equation (1)". 

2. Build staggered grid volumes, areas into geomechanical arrays — block 50b 

When the "input file" comprising the "mechanical properties of the rock" is entered, 
the staggered grid can be built. In this case, we need the staggered grid 'volumes' 
(how. much volume space the staggered grid would enclose) and the 'areas' (how 
much surface area the staggered grids would include). These 'volumes' and 'areas' 
are built into special 'arrays* used for geomechanics calculations. 

3. Build special NNC's for staggered srid in NNC stmcture — block 50c 

When the aforementioned 'arrays', used for the geomechanics calculations, are built, 
the next step is to build special NNC's. The acronym NNC refers to "Non-Neighbor 
Connections". Because the above defined 'equation (1)' [the "rock momentum 
balance equations"] are more difficult to solve than the standard reservoir simulation 
equations, it is necessary to create extra connections on the different grid blocks 
during the simulation in order to solve 'equation (1)' properly. Therefore, during this 
step, we build up these "Non-Neighbor Connections" for the staggered grid in a 
special "Non-Neighbor Connections Structure" or "NNC stmcture" that has been 
created inside the simulator. The "NNC stracture" is a collection of arrays and special 
areas that have been set aside inside the simulator. 
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4. Build boundary structure based on staggered gridding - block 50d 

When the ^^NNC structure" has been created inside the simulator, it is necessary to 
build up the boundaries representing the edges of the simulation or grid. The edges of 
the grid need this 'boundary structure', which is, in turn, based on the staggered grid. 
The whole *boundary structure*, which is a collection or arrays representing the 
'boundary conditions', has to be set up in order to rigorously solve the above 
mentioned equations, namely, the "rock momentum balance equations" defined above 
as 'equation (1)'. 

5. Assign input mechanical parameters to staggered grid array structure — block 50e 

When the 'boundary structure' has been build based on the staggered gridding, it is 
necessary to assign input mechanical parameters to the staggered grid array stmcture. 
After this is accomplished, the time stepping can now commence. 

6. Start time stepping - block 50f 

When the input mechanical parameters are assigned to the staggered grid array 
structure, the simulator will now start 'time stepping'. The simulator 46b of figure 15 
will time step forward in time and project into the future to determine what will 
happen. 

7. Calculate and store equation residuals for rock and fluid force balance in staggered 
grid - block 50g 

The first thing that is accomplished after the 'time stepping' is commenced is to 
calculate and store the equation residuals for the 'rock and fluid force balance'. The 
'rock and fluid force balance' refers to 'equation (1)', that is, the "rock momentum 
balance equations", which are set forth above and are duplicated below as follows: 
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The simulator wiD actually calculate the "rock momentum balance equations" 
representing 'equation (1)% and the residuals are exactiy as set forth in 'equation (1)'. 
Referring to 'equation (1)' set forth above (that is, the "rock momentum balance 
equations"), the term "residual" can be defined as follows: "How different from '0' is 
the left-hand side of the above referenced 'equation (1)'?". The objective is to try to 
make the left-hand side of 'equation (1)' equal to '0*, which is the right-hand side of 
the 'equation (1)\ The simulator 46b of figure 15 will therefore "calculate the 'rock 
momentum balance equations' representing 'equation (l)',.that is, the simulator 46b 
will adjust the parameters (u, v, w) and porosity and other resultant variables as set 
forth above in equations (2) and (3), until the left-hand side of the above 'equation 
(1)' is equal to the right-hand side of 'equation (1)', where the right-hand side of 
'equation (1)' is equal to *0'. Therefore, we know that we have solved the above 
referenced "rock momentum balance equations" of 'equation (1)' when the left-hand 
side of 'equation (1)' is equal to '0\ We find out how far away from '0' we are by 
calculating the 'residuals'. 



8. Calculate and store equation residual derivatives for rock and fluid force 
balance - block 50h 

At this point, we calculate and store the equation residual derivatives (which are 
required for the Newton gradient search), which derivatives will drive the "rock 
momentum balance equations" residuals [identified as 'equation (1)'] to '0\ At this 
point, we have: (1) 'residuals', and we have (2) 'derivatives of these residuals'. 
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9. Solve this linear system - block 50i ^ 
Now that we know: 

(1) the 'Residuals', and 

(2) the 'Derivatives of these Residuals', 

we can now solve this whole 'system of linear equations* representing the 'rock 
momentum balance equations* of 'equation (1)' which are set forth again below, as 
follows: 

ox oy az 

t^^*^-^-'- "> 

The above referenced 'rock momentum balance equations' of 'equation (1)' are 
solved together and simultaneously with the standard reservoir simulation equations. 

When the above referenced 'rock momentum balance equations' of 'equation (1)' are 
solved, the "u", "v", and "w" displacement parameters, representing movement of the 
rock, are determined, where: 

u = rock displacement in the x-direction 
V = rock displacement in the y-direction 
w = rock displacement in the z-direction 

When the displacement parameters (u, v, w) are determined, the above referenced 
'equation (3)' is used to determine (the 'x,y,z elongation strains') and 
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"^xyo'z.zx" (the *shear strains'), since "£,,3,,," and are function of "u", "v", 

and "w", as follows: 



dx 



''-By 
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When "£,yj" and "7xy.yz,zx" ^® determined from 'equation (3)% the above 
referenced 'equation (2)' is used to determine "o,.y.j " (the 'elastic normal rock stress 
in x,y,z directions') and "i^y^^" (the 'elastic shear stress') since "0^.^^" and 
"■^xy-yz-xt " are a function of " fi, " and " " in 'equation (2)' as follows: 



o,=2G£, + 4£x + £,+«J 



(2) 



When the "o^^^" and "ixy,yz.xz" ^® determined from 'equation (2)', the above 
referenced 'equation (1)' is solved, since 'equation (1)' is a function of "Oj^^^" and 
"-r^ as follows: 
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dx dy dz 




+ Fy + Py =0 



(1) 



The result of these calculations will produce the Table of figure 16. A reservoir 
engineer would be interested in knowing each of the quantities shown in the Table of 
figure 16. 

10. Evaluate Failure Criteria and Solve for Plastic Displacements — block 50j 

A final step is required if the reservoir engineer has specified regions in the step 50a 
in Figure 18 where the rock may fail. In this case, the final "Oj, ^.^ " and "t^^,^^^'* are 

used to evaluate the chosen failure criterion which may be Mohr-Coulomb as referred 
to above in the section 'Description of the Invention' subsection Plasticity. If the 
failure criterion is exceeded, then a set of equations is set up to solve for the plastic 
displacements, denoted (Up, Vp, Wp) above, which will force the failure criterion 
residual to be exactly zero. This system is set up and solved exactly as described 
above with the exception that the residuals are now the failure criterion residuals 
which are driven to zero instead of the equation (1) residuals. This step is 50j in 
Figure 18. 

In figure 18, when the *system of linear equations' set forth above as 'equation (1)' 
has been solved, referring to the 'feed-back loop' 64 in figure 18, the simulator 46b of 
figure 15 continues its 'time stepping* by repeating the implementation of blocks 50g, 
50h, 50i and 50j in figure 18. 

When the flow chart of figure 18 is executed, the 'result* is: a solution to a problem, 
which a reservoir engineer would want to solve, in which you have certain boundary 
conditions applied to the rock, boundary conditions applied to 'equation (1)', and 
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boundary conditions applied to the reservoir simulation equations. For example, in 
the North Sea, problems with 'subsidence' exist. The term 'subsidence' is a situation 
where the ground actually 'gives-way' or 'fails'. The above calculation would be 
applied to this 'subsidence' problem. The above type of calculation, which represents 
a solution to the above referenced 'Elastic Stress Equations' using 'staggered 
gridding', would be used to handle the 'more difficult' fields, such as the fields having 
the 'subsidence' problem of fiigure 17. 

In summary, rock will initially behave elastically because it is in equilibrium and there 
are no extemal forces. Elastic behavior means that the rock behaves like a spring; that 
is, you can push it (i.e., apply force) and it will return exactly to its original position 
when you stop pushing. Subsidence happens when the rock is pushed beyond its own 
strength. For example, in a particular reservoir where certain wells are withdrawing 
fluids, the rock may no longer have the support of the water pressure, and the weight 
of the overburden will become high enough to cause the rock to fail. If the rock fails, 
e.g. by crushing or cracking, then, it no longer moves elastically but rather moves 
plastically. This plastic movement is what is referred to as "subsidence". We know 
this will occur because we solve 'equation (1)', which are the rock momentum 
balance equations (otherwise known as rock force balance equations). If these 
equations tell us that the forces become too large in one direction relative to the forces 
in Another direction at some particular point in time, then, we know that failure will 
occur. At this point, we further solve a final set of equations which will calculate the 
plastic displacements. "Derivatives" are needed because these equations are 
nonlinear. 

References 
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The invention being thus described, it will be obvious that the same may be varied in 
many ways. Such variations are not to be regarded as a departure from the spirit and 
scope of the invention, and all such modifications as would be obvious to one skilled 
in the art are intended to be included within the scope of the following claims. 
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CLAIMS 

1. A method of determining subsidence in a reservoir, said subsidence resulting from 
a continual withdrawal of hydrocarbon deposits and other fluids from said reservoir, 
comprising the steps of: 

(a) estimating rock displacement parameters u, v, and w representing rock movement 
in the x, y, and z directions; 

(b) determining e^ ^^ (%y,z elongation strains') and J^y^y^^zx ('shear strains') from the 
rock displacement parameters u, v, and w; 

(c) determining o^ ^ ^ ('elastic normal rock stress in x,y,z directions') and t^ y^ ^^ 
('elastic shear stress') from the s^ y ^ and J^^^^y^^^ ; 

(d) solving a set of rock momentum balance differential equations from the o^ y ^ and 
the 'txy.yz.xz detemaining if any residuals exist, the rock displacement parameters u, 
V, and w representing accurate rock displacement parameters for said reservoir when 
the residuals are substantially equal to zero; and 

(e) determining said subsidence in said reservoir from said accurate rock displacement 
parameters determined from step (d), by solving a final set of failure criterion 
equations, which comprise the residuals and any derivatives, to determine a set of 
rock plastic displacements, said rock plastic displacements forming a part of said rock 
displacement parameters u, v, and w. 
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2. The method of claim 1, wherein the detemiining step (b) comprises the step of: 

deteraiining said £^ (the *x,y,z elongation strains') and said "i^^y^^^ (the 'shear 

strains') from said rock displacement parameters u, v, and w by using the following 
equations: 
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3. The method of claim 2, wherein the determining step (c) comprises the step of: 

determining said o^ y ^ (the 'elastic normal rock stress in x,y,z directions') and said 
'elastic shear stress') from said and 7xy.yz.zx using the following 

equations: 
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4. The method of claim 3, wherein the solving step (d) comprises the steps of: 

solving said set of rock momentum balance differential equations from said o^ ^ ^ and 
said t 

jcy.yz.xz » 



wherein said rock momentum balance differential equations have a left-hand side on a 
left side of an equal sign and a right-hand side on a right side of said equal sign, said 
rock momentum balance differential equations including. 



^ + -:r^ + ^^ + F.+P. =0 



dx dy dz 
dcT^ 3t^. dr^ 



dy dx dz 



da^ dr^ dr. 



dz dx dy ' ' 



said rock displacement parameters u, v, and w estimated in step (a) representing said 
accurate rock displacement parameters for said reservoir when a left-hand side of said 
rock momentum balance differential equations substantially equal a right hand side of 
said rock momentum balance differential equations, 

5. The method of claim 4, wherein the solving step (e) comprises the steps of using 
the said o^^y^ and said "txy^y^^xz determine whether a failure criterion has been 
exceeded in one of a set of staggered grids. 
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6. The method of claim 5, wherein, when said failure criterion has been exceeded in 
said one of said set of staggered grids, the solving step (e) further comprises the step 
of: 

setting up a final set of equations, each equation being a failure criterion 
approximately equal to zero, and 

determining a set of parameters up, vp, wp which are plastic displacements resulting 
from a failure of the rock, said plastic displacements forming a part of said subsidence 
and being determined from said accurate rock displacement parameters for said 
reservoir. 

7. A program storage device readable by a machine, tangibly embodying a program 
of instructions executable by the machine to perform method steps for deternodning 
subsidence in a reservoir, said subsidence resulting from a continual withdrawal of 
hydrocarbon deposits and other fluids from said reservoir, said method steps 
comprising: 

(a) estimating rock displacement parameters u, v, and w representing rock movement 
in the x, y, and z directions; 

(b) determining £^ ^ ^ Cx»y»z elongation strains') and y^y^y^^^ ('shear strains') from the 
rock displacement parameters u, v, and w; 

(c) determining o^ y ^ (*elastic normal rock stress in x,y,z directions') and t^ ^^^ 
('elastic shear stress') from the e^ y ^ and }xy,yi,zK » 

(d) solving a set of rock momentum balance differential equations from the o^ y ^ and 
^xy,yz,xz determining if any residuals exist, the rock displacement parameters u, 
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V, and w representing accurate rock displacement parameters for said reservoir when 
the residuals are substantially equal to zero; and 

(e) determining said subsidence in said reservoir from said accurate rock displacement 
parameters determined from step (d), by solving a final set of failure criterion 
equations, which comprise the residuals and any derivatives, to determine a set of 
rock plastic displacements, said rock plastic displacements forming a part of said rock 
displacement parameters u, v, and w. 

8. The program storage device of claim 7, wherein the detennining step (b) 
comprises the step of: 

determining said e^ y^^ (the 'x,y,z elongation strains') and said J^^y^^^ (the 'shear 

strains') from said rock displacement parameters u, v, and w by using the following 
equations: 





_ 3m 






dx 






dv 












dw 




" ai" 






du 






dy 






dv 


dw 


dz 


dy 




_ dw 


du 
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9. The program storage device of claim 8, wherein the determining step (c) comprises 
the step of: 

determining said o^ y^^ (the 'elastic normal rock stress in x,y,z directions') and said 
"txy.yz^xz *elastic shear stress') from said and J^^^y^^^ by using the following 

equations: 



(Ty^lGey + Xif^ + Ey + e^) 

(73=2G£,+4^, + £3,+£j 



10. The program storage device of claim 9, wherein the solving step (d) comprises 
the steps of: 

solving said set of rock momentum balance differential equations from said o^ y ^ and 
said 

wherein said rock momentum balance differential equations have a left-hand side on a 
left side of an equal sign and a right-hand side on a right side of said equal sign, said 
rock momentum balance differential equations including. 



— ^ + F, =0 

ax 83; 9z ' ' 

da^ dr^ dr^ 
— ^-h — ^-h — ^ + F +P =0 

dy dx dz ' ' 

da^ dr^ „ ^ ^ 

dz dx dy ^^-^ ^ 
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said rock displacement parameters u, v, and w estimated in step (a) representing said 
accurate rock displacement parameters for said reservoir when a left-hand side of said 
rock momentum balance differential equations substantially equal a right hand side of 
said rock momentum balance differential equations, 

11. The program storage device of claim 10, wherein the solving step (e) comprises 
the steps of using the said o^ y ^ and said "t^y^yz^xz determine whether a failure 
criterion has been exceeded in one of a set of staggered grids. 

12. The program storage device of claim 11, wherein, when said failure criterion has 
been exceeded in said one of said set of staggered grids, the solving step (e) further 
comprises the step of: 

setting up a final set of equations, each equation being a failure criterion 
approximately equal to zero, and 

determining a set of parameters up, vp, wp which are plastic displacements resulting 
from a failure of the rock, said plastic displacements fomiing a part of said subsidence 
and being determined from said accurate rock displacement parameters for said 
reservoir. 
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